function answer = vlasov(plist)
   pcopy = plist;
   Tx = plist(3);
   n =  plist(6);
   C =  plist(7); %%C*v_temp = v_max
   v_temp = 4.19e7*(Tx^0.5);
   v_max = C*v_temp;
   dv = v_max/n;
   v_arr = -v_max:dv:v_max;
   f0 = create_f0(plist);
   t_start = plist(8);
   t_end   = plist(9);
   B0 = plist(1);
   e = 4.8032043e-10;
   m = 9.1093822e-28;
   c = 3e10;
   w_ce = e*B0/(m*c);
   dt = 0.01*2*pi/w_ce;
   pcopy(11)=dt;
   for t = t_start:dt:t_end
       pcopy(10) = t;
       T = operator_matrix(pcopy);
       f1 = f0\T;
       f1_mat = reshape(f1, 2*n+1, 2*n+1);
       [vx,vy]=meshgrid(v_arr, v_arr);
       p=pcolor(vx, vy, f1_mat);
       set(p,'LineStyle', 'none');
       f0 = f1';
       pause(100e-3);
   end
   answer = 0;
   
    
function [M, N] = create_block(plist, nblock, A, B)
    n = plist(6);
    nu = plist(13);
    dt = plist(11);
    size = 2*n+1;
    M = sparse(size, size);
    N = sparse(size, size);
    N(1,1) = 0;
    N(end, end) = 0;
    M(1,1) = 1;
    M(end,end) = 1;
    for i = 2:1:(size-1)
        %%%This block can be faster!
        M(i, i) = (1+A(i) + B(nblock))/(1-nu*dt);
        M(i, i + 1) = -B(nblock);
        N(i, i) = -A(i);
    end
    
function T = operator_matrix(plist)
    e = 4.8032043e-10;
    m = 9.1093822e-28;
    c = 3e10;
    B0 = plist(1);
    E0 = plist(2);
    Tx = plist(3);
    %Ty = plist(4);
    T =  plist(5);
    n =  plist(6); %%(2*n+1)X(2*n+1) - is size of operator's approx. matrix
    C =  plist(7); %%C*v_temp = v_max
    %t_start = plist(8);
    %t_end   = plist(9);
    t_curr  = plist(10);
    dt = plist(11);
    %n0 = plist(12);
    
    
    %%suppose Tx == Ty
    v_temp = 4.19e7*(Tx^0.5);
    v_max = C*v_temp;
    dv = v_max/n;
    w_ce = e*B0/(m*c);
    D = dt/dv;
    e_m = e/m;
    E = E0*exp(-10*(t_curr+dt)*w_ce)*sin(w_ce*(t_curr+dt))*e_m;
    
    v_arr = -v_max:dv:v_max;
    size = 2*n+1;
    T = sparse(0,0);
    m_row = sparse(0,0);
    I_sparse = sparse(eye(size,size));
    O_sparse = sparse(zeros(size,size));
    for j = 2:1:size
        m_row = horzcat(m_row, O_sparse);
    end
    T = horzcat(I_sparse, m_row);
    B = zeros(size,1);
    A = zeros(size,1);
    
    B = -D*w_ce*v_arr;
    A = D*(w_ce*v_arr - E);
    m_row = sparse(0,0);
    for i = 2:1:(size-1)
        
        [M, N] = create_block(plist, i, A, B);
        m_row = sparse(0,0);
        for j = 1:1:size
            
            if (j < i)
                m_row = horzcat(m_row, O_sparse);
            else
                if (i == j)
                    m_row = horzcat(m_row, M);
                else
                    if (j == i + 1)
                        m_row = horzcat(m_row, N);
                    else
                        if (j > i + 1)
                            m_row = horzcat(m_row, O_sparse);
                        end
                    end
                end
            end 
        end
        T = vertcat(T, m_row);
    end
    m_row = sparse(0,0);
    for j = 2:1:size
        m_row = horzcat(m_row, O_sparse);
    end
    m_row = horzcat(m_row, I_sparse);
    T = vertcat(T, m_row);
    
 function f0 = create_f0(plist)
    Tx = plist(3);
    Ty = plist(4);
    n =  plist(6);
    C =  plist(7); %%C*v_temp = v_max
    v_temp = 4.19e7*(Tx^0.5);
    v_max = C*v_temp;
    dv = v_max/n;
    v_arr = -v_max:dv:v_max;
    n0 = plist(12);
    m = 9.1093822e-28;
    
    A = m/(2*1.6e-12*Tx);
    size = (2*n+1);
    f0 = zeros(size^2,1);
    
    for i = 2:1:(size-1)
        for j = 2:1:(size-1)
            f0((i-1)*size+j) = n0*(A/pi)*exp(-A*(v_arr(i)^2+v_arr(j)^2));
        end
    end
    
        
        
